Critical Correlations for Short- Range Valence-Bond Wave Functions on the Square Lattice 
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We investigate the arguably simplest SU(2) -invariant wave functions capable of accounting for spin-liquid 
behavior, expressed in terms of nearest-neighbor valence-bond states on the square lattice and characterized by 
different topological invariants. While such wave-functions are known to exhibit short-range spin correlations, 
we perform Monte Carlo simulations and show that four-point correlations decay algebraically with an exponent 
1.16(4). This is reminiscent of the classical dimer problem, albeit with a slower decay. Furthermore, these 
correlators are found to be spatially modulated according to a wave- vector related to the topological invariants. 
We conclude that a recently proposed spin Hamiltonian that stabilizes the here considered wave-function(s) as 
its (degenerate) ground-state(s) should exhibit gapped spin and gapless non-magnetic excitations. 

PACS numbers: 75.10.Kt,75.10Jm,75.40.Mg 



Introduction — The quest for quantum spin-liquid (QSL) 
states of matter 1 is a longstanding research topic that can 
be traced back to Anderson's proposal. 2 Building on earlier 
work, he conjectured that strong quantum fluctuations, en- 
hanced by frustration and/or low coordination, would weaken 
S'?7(2)-broken order and occasionally drive an antiferromag- 
net towards a "disordered" state with exponentially decaying 
spin correlations, describable in terms of short-ranged spin- 
singlet, or valence-bond (VB), degrees of freedom. 2 3 Interest 
in Anderson's insight was further triggered in connection with 
the cuprates since, soon after their discovery, spin-singlets in 
a QSL were interpreted as "pre-formed Cooper pairs" that 
would superconduct upon doping. 4 

Major advances (reviewed in Refs. 1, 5, and 6) have taken 
place since the original proposal, 2 including a classification of 
possible QSL states 7 and explicit realizations in lattice mod- 
els in dimension d > 1. Also, theoretical ideas put for- 
ward in the early days of high-T c 1718 have been considerably 
developed and resulted in a full-fledged formalism 19 as well 
as efficient numerical approaches 20 for handling VB states. 
On the experimental side, a number of compounds have been 
shown not to display magnetic order down to the lowest ac- 
cessible temperatures, 1 much below the energy scale set by 
exchange interactions, and are thus candidates for the realiza- 
tion of QSLs. However, in spite of such advances, a complete 
characterization of QSL states is still missing, precluding un- 
ambiguous identification of experimental realizations, since 
absence of magnetic order does not exclude the occurrence of, 
for instance, more conventional valence-bond crystals (VBC) 
that break lattice symmetries (see e.g. Ref. 21). 

Within this context, we investigate a family of nearest- 
neighbor VB (NN-VB) states on the square lattice by per- 
forming Monte Carlo (MC) simulations based on a recently 
introduced algorithm. 22 We revisit the pioneering work by 
Sutherland, 18 where a closely related NN-VB state was in- 
troduced, and provide a thorough characterization of the ar- 
guably simplest S'J7(2)-invariant wave functions capable of 
accounting for QSL behavior. Although it has long been 
known that NN-VB states on the square lattice are non- 
magnetic, 17 the possibility of other types of order, such as 
VBC, has not yet been excluded. Despite their simplicity, 
and consequent theoretical appeal, the NN-VB states possess 



highly non-trivial properties, that we explore in what follows. 

Wave functions — NN-VB states are obtained by contract- 
ing spins attached to NN sites i and j of a lattice into a singlet 
state, [i,j] = ^(\t%ij) - \Utj))- Since each s P in onl y 
pairs with one of its neighbors at a time, there is a one-to-one 
correspondence between NN-VB configurations and those of 
hard-core classical dimers on the same lattice. 23-26 

A crucial property of VB states is their non-orthogonality. 
The overlap between two VB configurations is given by 
(V>i|V>2) = ±2 W£ -f , 18 where N = L 2 is the number of 
sites and Nc the number of loops in the transition graph ob- 
tained by superposing the dimer configurations associated to 
\tf>i) and \ip2) [Fig. 1(a)]. For the square and other bipartite 
lattices, that can be split into two sublattices A and B, overlaps 
between arbitrary VB states can be ensured to be always posi- 
tive, so that stochastic methods apply (see below), by choosing 
i G A and j G B in the anti-symmetric singlet 

An additional important point concerns the fact that dimer 
configurations can be split into topological sectors. 6 On a 
torus, the transition graph for two dimer coverings belonging 
to different topological sectors displays non-local loops that 
wind around the system [Fig. 1(a)], so that one dimer config- 
uration can not be continuously deformed onto the other via 
local dimer rearrangements. For bipartite lattices the num- 
ber of topological sectors is extensive and each sector can be 
labelled by topological invariants termed winding numbers, 
w = (w x , w y ): w x (w y ) is defined as the difference between 
the number of B <— A and A — > B dimers along a reference 
line in the y (x) direction [Fig. 1(a)]. VB configurations char- 
acterized by different w are orthogonal in the thermodynamic 
limit: the transition graph between two such configurations, 
c Wl ) and |c W2 ), contains at least one winding loop that com- 
prises a minimum of L dimers, so that Nc < (N — L)/2, im- 



plying that 



|c W2 ) < 2 2 and vanishes when L — > oo. 



Having introduced the ingredients, we are able to write 
down the NN-VB wave functions we wish to investigate: 



(i) 



In contrast to the wave function analyzed by Sutherland, 18 
who did not take the existence of topological sectors into ac- 



FIG. 1. (Color online) (a) Transition graph between two NN-VB 
configurations on the square lattice (sublattices AI3 are indicated by 
filled/open circles). Reference lines for the winding numbers w = 
(w x ,w y ) are indicated by dashed lines: configurations with w = 
(1,0) (red lines) and w = (0, 1) (blue lines) are shown. Trivial 
loops with coinciding VBs are depicted as thick-black lines and a 
loop winding in both directions is evident, (b) NN spin correlations 
versus L^ 1 (from MC simulations) for winding sectors w = (0, 0), 
(0, 1), (1, 0) and (1, 1). For w (0, 0), correlations along ±e x and 
±e y are discriminated. Lines are linear-quadratic fits. 

count and considered an equal amplitude superposition of all 
NN-VB states, each \t/j w } is an equal amplitude superposi- 
tion of VB configurations |c w ) with fixed winding numbers 
w = (w x , w y ). Our motivation for doing so is our previous 
remark that (VvjVVa) = f° r w i 7^ w 2 in the thermo- 
dynamic limit. Interestingly, this implies that each |i/) w ) is a 
(degenerate) ground-state on a torus of the local spin Hamil- 
tonian recently proposed by Cano and Fendley. 28 Although 
the number of winding sectors is extensive on a torus, we will 
show that it is possible to infer the properties of arbitrary states 
IV'w) by focusing on a few sectors with low w. 

Algorithms — The expectation value of an observable O in 
Eq. (1) can be measured by evaluating 

itn\ 1 (c w ,i|0|c w 2 ) , , . 
C w = -= > —, ; r C w ,i Cw,2 , (2) 

Z ^ (Cw.i Cw, 2 ) 

C w ,i,C w .2 1 ' ' ' ' 

where Z = (VAvI^w) is the normalization. As first pointed 
out in Ref. 17, (C) w can be efficiently computed in a stochas- 
tic manner: the estimator (c Wi i|0|c Wj 2)/(c w ,i|c w ,2), that 
for most observables of interest is readily evaluated by an- 
alyzing the loop structure in the transition graphs, 19 is sam- 
pled by generating pairs of NN-VB configurations |c w .i) and 
|c w ,2) according to the statistical weight given by their overlap 
(c w ,i|cw,2) = 2 Ar ^ w;1 < 2 )-f [iV £ (w; 1,2) denotes the num- 
ber of loops in the transition graph (c w i|c Wj 2}]. 

Major advances in sampling techniques have been achieved 
since the work by Liang et al. 17 and particularly well suited to 
our purposes is a recently introduced algorithm. 22 Basically 
(we refer to Ref. 22 for details), one combines non-local up- 
dates for the underlying dimer configurations, so to ensure 
small auto-correlations times, with spin updates that allow 
for an efficient sampling of the overlap weight, with unitary 
acceptance rate. By relying on this algorithm, we simulate 
systems with periodic boundary conditions (PBC) of sizes 
of up to L = 128. Topological symmetry is easily imple- 
mented in the simulations by starting from a configuration in 



FIG. 2. (Color online) (a) Spin correlation ( — l) r (S ■ S r ) between 
the spin at the origin and spins located at r = xe x + ye y . (b) 
(— l) r (So ■ S r ) versus distance r. An exponential fit, ( — l) r (So • 
S r ) ~ exp(— r/£), yields the correlation length £ = 1.35(1). Data 
fori = 128 and w = (0,0). 



a given winding sector and discarding measurements for all 
MC moves that change w: for large L winding updates are 
exponentially rare and this only causes small efficiency losses. 

Short-range spin order — We start by analyzing the spin 
texture in wavefunction Eq. (1). NN spin correlations, (S r • 
S r ±e a ) (e a is the unit vector in the a = x,y direction) are 
plotted as a function of inverse system size L _1 in Fig. 1(b) 
for w = (0,0), (1,0), (0,1) and (1,1). Deviations among 
results for different w are observed for small systems and, 
additionally, for w ^ (0, 0) vertical/horizontal and A — >• B 
I B — >• A correlations differ. However, all data converge to 
the same value in the thermodynamic limit, according to a 
linear-quadratic best-fit analysis. Although we have no rigor- 
ous justification for such scaling, we obtain (S r • S r + e „ ) = 
—0.295953(7) when L — >• oo for all sectors. Spin correla- 
tions (So • S r ) as a function of distance are plotted in Fig. 2(a- 
b) for L = 128 and w = (0, 0) (virtually identical results 
are obtained for other L and w) and display a perfect stag- 
gered pattern consistent with (isotropic) short-range Neel or- 
der. (So • S r ) decays very fast with |r| and an exponential fit 
[Fig. 2(b)] yields £ = 1.35(1) for the spin correlation length. 

Critical correlations — We proceed to the characterization 
of "dimer order" by analyzing the four-point connected cor- 
relators C^ kl = ((Si • 8,0(8* • SO) - (Si ■ Sj)(S k • S ; ), 
where both i,j and k, I are NN sites on the square lattice. In 
Fig. 3 we show MC data for the spatial dependence of rC l ^ u 
(r is the distance between dimers) for L — 16 and sectors 
w = (0, 0), (1, 0), (0, 1) and (1, 1). We first notice that both 
Cfp (correlations for parallel dimers i, j and k, I) and C^ 1 
(perpendicular dimers i, j and k, I) are spatially modulated 
for w 7^ (0, 0). Inspection of the results in Fig. 3, and similar 
ones for higher w (not shown), allows us to deduce that mod- 
ulation for Cm 7 ' [C 1 ^ 1 ] is entirely accounted for by a phase 
factor cos(Q-r) [sin(Q-r)], with a wave-vector given in terms 
of the winding numbers, Q = =£-{w v , w x ). This inference is 
confirmed by our quantitative analysis below. 30 Furthermore, 
we notice that no clear spatial dependence is noticeable for 
rC l ^ kl in Fig. 3(a), suggesting that four-point correlations in 
Eq. (1) decay algebraically with r with an exponent close to 
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(c)-w = (0, 1) (d)-w = (l, 1) 




FIG. 3. (Color online) MC results for rC ljkl , for L = 16 (PBC) and 
indicated w. In all panels, the reference bond is indicated by a thick- 
black line and the thickness of the remaining ones is proportional 
to rC y : blue (pale-red) lines indicate positive (negative) values. 
Lines along which CtP is strongest are indicated by black circles. 
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FIG. 4. (Color online) MC results for Cfi kl as a function of distance, 
(a) Longitudinal correlations C^ kl for w = (0, 0) and various sys- 
tem sizes, (b) C^ kl along zero-phase anti-nodal directions (filled cir- 
cles in Fig. 3) for L = 128 and w = (0, 0), (1, 0), (0, 1) and (1, 1). 
(c) C^ kl (r) I cos(Q • r) and C^' (r) / sin(Q • r) for w = (0, 1) and 
L — 128 (data separated from a nodal line by a parallel displacement 
dx < 8 are excluded), (d) Absolute value of transverse correlations 
C*f for w = (0, 0) and L = 128, along the line highlighted in 
Fig. 3(a). In (a-c) the line indicates our best fit yielding the exponent 
a = 1.16(4) and in (d) the fit yielding a' = 2.53(5). 



unity (see below). 

In Fig. 4(a-b) we plot C^ kl (r) along the anti-nodal lines 
with strongest correlations (hence smallest relative errors), 
highlighted in Fig. 3, for: (a) w = (0,0) and various sys- 
tem sizes L and (b) L = 128 and w = (0,0), (1,0), 
(0, 1) and (1, 1). Results are consistent with power-law de- 
cay, Cm 3 ' ~ r~ a , as conjectured in Refs. 18 and 31. Fitting 
the data in Fig. 4(b) we arrive at a = 1.16(4). Although small 
deviations from algebraic behavior are seen for large distances 
in Fig. 4(a-c), the value of r at which they start to occur in- 
creases linearly with L (not shown) and we thus conclude that 
this "upturn" in Fig. 4(a-c) is merely a finite-size effect. 

We analyze the spatial modulations for correlations and in 
Fig. 4(c) we plot &J k \r) and Cf kl {r) for L = 128 and 

w = (0,1). By respectively dividing OfP (r), C l [ kl (r) 
by the phase factors cos(Q x dx), sm(Q x dx) [dx = along 
the anti-nodal line for C 1 ^ (r) highlighted in Fig. 3(c) and 
Q = ^(1, 0) in this case], we notice that all curves collapse, 
confirming that C$ (r) and C l / kl (r) decay with the same 
exponent and are indeed modulated according to the phase 
factors cos(Q • r) and sin(Q • r), with Q = ^(w y ,w x ). 

Sub-leading corrections to the scaling exponent can be ob- 
tained by analyzing (ff in the w = (0, 0) sector [see 
Fig. 3(a)]. Data for (ff are plotted as a function of r 
in Fig. 4(d), for L = 128 and w = (0,0). A fit yields 
a' = 2.53(5) for the sub-leading exponent. 

Finally, we address the point of what specific type of quasi- 
long-range dimer order is encoded in Eq. (1). In doing so, 
we analyze the dimer order parameter D defined by D a = 
N" 1 X) r (~ l) r °Sr ■ S r+ecl . Due to the absence of long-range 
order, (D) is expected to vanish when L —> oo. However, 
information concerning the symmetry of the quasi-ordered 
state is obtainable by analyzing the angular dependence on 
4> = arctan(Dj / /D a; ) in the histogram P(D x ,D y ) for oc- 
currences of D x and D y in the simulations. In Fig. 5 we 
plot P(D x ,Dy) for L = 96 and w = (0,0). Commonly 
observed VBCs on the square lattice, "columnar" and "pla- 
quette" states (for a detailed account see Ref. 21), would dis- 
play, respectively, peaks located at <f> = {0, ±7r/2,7r} and 
4> = {±7r/4, ±37r/4}. However, no angular structure is ev- 
ident in Fig. 5 and data are in favor of an continuous U(l) 
symmetry, a priori different from the U(l) symmetry associ- 
ated to topological degeneracy. 

Conclusions — We have investigated NN-VB wave func- 
tions on the square lattice [Eq. (1)], characterized by wind- 
ing numbers w, by performing MC simulations. We con- 
firm earlier findings in favor of short-ranged spin order 17 ' 29 
and, more interestingly, we find that dimer-dimer correla- 
tions are critical, a situation reminiscent of the one encoun- 
tered for classical dimers 26 and thus for the ground-state of 
the quantum dimer model (QDM) on the square lattice at the 
Rokhsar-Kivelson point. 32 However, such correlations decay 
considerably slower for Eq. (1) than in the classical case, sug- 
gesting increased tendency towards VBC order: an exponent 
a = 1.16(4) accounts for the decay of both longitudinal and 
transverse correlation for all w studied, while in the classical 
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FIG. 5. (Color online) Histogram P(D X , D y ) for L = 96 and w = 
(0,0). 

case one has a c i ass = 2. 26 In this context, it would be in- 
teresting to analyze how exponents evolve by considering the 
overlap as a tunable parameter, 32 between the herein studied 
case and the limit of orthogonal dimer configurations of the 
QDM. 32 

From a broader perspective, we analyze how the wave- 
function Eq. (1) fits into the general classification of QSL 
states. 7 While we stress that our work concerns wave- 
functions and not the full spectrum of a given Hamiltonian, 
we notice that each |^) w ) [Eq. (1)] is a (degenerate) ground- 
state of the local model of Ref. 28 on a torus. We thus predict 
such 5?7(2)-invariant Hamiltonian to display gapped spin ex- 
citations, due to short-ranged spin order, 18 ' 33 and gapless non- 
magnetic excitations, since four-point correlations are criti- 



cal and a theorem by Hastings applies. Adopting the termi- 
nology of Ref. 7, the latter excitations correspond to gapless 
gauge modes. Given its extensive degeneracy on a torus, our 
results altogether suggest that the ground-state of the model 
of Ref. 28 is a gapped U(l) or SU(2) spin liquid, 7 a state be- 
lieved to be unstable for a generic local spin model. In this 
context, the complete characterization of the spectrum of the 
model of Ref. 28 and of perturbations thereof would be of high 
interest. 

Directions for further research opened up by our work in- 
clude the study of wave-functions similar to Eq. (1) in dif- 
ferent geometries. In d = 3, we expect NN-VB states on 
bipartite lattices to display spin order, as it happens for the 
simple cubic lattice, 29 but it would be interesting to search 
for traces of the Coulomb phase of dimer models in d = 3. 35 
NN-VB wave-functions on d = 2 geometrically frustrated lat- 
tices also deserve investigation. 1213 However, a sign problem 
precludes MC simulations as we perform here and alternative 
approaches are called for in this case. 

Note Added — While preparing this manuscript we became 
aware of related work by Tang et al.. 36 
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